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I .  '^x'ly^ls  of  a  Shutdcum  Queueing  Svste'-j. 

1.  Introduction.  ^  The  pumose  of  this  study  is  to  analyze  a  special  type 
of  job  shop  queueing  system  which  has  the  following  features. 

(i)  A  finite  number  of  customers  is  present  initially  and  no  new 
customers  arrive; 

(II)  there  are  two  service  stations,  each  with  Its  oi^n  Input  and 
each  serving  the  output  of  the  other;  and 

(III)  as  soon  as  a  customer  is  served  by  both  stations,  he  leaves 
the  system. 

A  queueing  system  with  these  features  is  called  a  shutdown  queueing  system 
and  it  may  also  be  called  a  clearance  problem. 

Queueing  problems  of  this  type  arise  in  connection  with  fuel  and  supply 
facilities  for  aircraft.  A  fixed  number  of  aircraft  are  scheduled  for  fueling 
end  supply  loading.  .  "  .  -  /  ,  ■  , 

Each  of  these  operations  is  conducted  at  two  separate  stations  called 
station  I  and  station  II.  Initially,  station  I  (which  may  be  the  fuel  station) 
has  M  customers  while  station  II  (supply  station)  has  N  customers.  After  an 
aircraft,  henceforth  called  a  customer,  initially  at  station  I  completes  its 
service,  it  enters  the  queue  for  service  at  station  II.  When  Its  service  at 
station  II  is  also  complete,  it  leaves  the  system.  The  service  time  at  both 
stations  is  assumed  to  be  a  random  quantity.  Slmilariy,  a  customer  Initially 
in  the  station  II  queue,  moves  to  the  station  I  queue  and  leaves  the  system 
after  his  service  Is  complete. 

The  length  of  time  required  for  the  entire  operation,  or  the  clearance 
of  the  queue,  Is  a  random  quantity  whose  probability  distribution  is  of  major 
interest.  Also  of  Interest  is  the  idle  time  distribution  at  each  of  the  supply 
stations.  Since  no  new  customers  join  the  queue,  only  nonsteady  state  results 
are  of  Interest. 

Under  the  assumption  .that  the  service  time  distrlbutipn  Is  negative 
exponential,  this  problem  was  studied  by  Milch  and  Waggoner 

They  obtained  Laplace  transforms  of  the  probability  distribution  of  the  total 
cneratlonal^lme,  the  idlS  time  distribution  for  each  station  and  the  waiting 
time  distribution  for  various  customers  In  the  system.  These  transforms  were 
quite  complex  and  not  easily  Inverted. 
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Thtj  work  of  llilch  o.nd  ’/aggoncr  can  be  studied  ;*.n  greater  detail  eed 
generality  using  discrete  time  methods  rather  than  continuous  time  models. 
C;’scretizing  the  time  variable  makes  it  possible  to  retain  all  conditional 
probability  distributions  for  the  entire  time  interval  of  the  operation  of  both 
stations.  Because  of  the  many  possible  interactions  between  the  two  queues,  the 
analysis  of  this  general  model  by  classical  analytical  methods  is  an  Intractable 
problem.  The  analysis  of  Mich  and  Waggoner  relies  heavily  on  the  memoryless 
property  of  the  exponential  41strlbutlon.i 

Even  using  discrete  time,  the  numberlcal  solution  of  this  problem 
requires  extensive  computing.  Many  realizations  of  the  process  occur  with 
^‘-gllgible  probability,  which  may  be  conveniently  discarded  early  in  the  complex 
computational  algorithm  reported  on  here.  On  the  other  hand,  it  is  possible 
to  obtain  upper  bounds  on  the  probability  neglected  at  each  stage. 

2.  Statement  of  the  Problem.  Consider  two  stations,  I  and  II.  Initially, 

station  I  has  N  customers  and  station  II  has  N  customers.  The  service  time  at 
each  station  is  measured  in  discrete  time  units.  The  probability  that  a 
customer  at  station  I  requires  i  units  of  service  time  is  for 

1  “  1,  2,  . . . ,  Lj^  while  at  station  II  a  customer  requires  j  units  of  service 
with  probability  Yj  for  J  ■  1,  2,  ...,  L2.  The  two  stations  operate 
independently.  After  a  customer  conq>lete8  service  at  the  station  of  his 
initial  assignment,  he  Joins  the  end  of  the  queue  of  the  other  station.  The 
amount  of  transit  time  la  assumed  to  be  zero.  (This  assumption  is  not 
restrictive,  since  the  transit  time  could  be  counted  as  part  of  the  service  time.) 

After  being  served  by  both  stations,  the  customer  leaves  the  system. 

It  is  possible  for  the  queue  at  either  station  to  become  empty.  The  station 
cr.en  becomes  Idle  and  remains  so  until  a  customer  completes  service  at  the 
other  statloh.  Such  Idleness  may  repeatedly  occur  during  the  service  operation 
of  the  system. 

3.  Random  Walk  Representation.  The  progress  of  this  queueing  system  may 
be  represented  by  a  random  walk  process  on  a  two  dimensional  lattice.  The 
random  walk  starts  at  time  zero  at  the  origin.  Thereafter,  the  system  is 


Figure  1  Random  Walk  Representation 


represented  by  the  pair,  (v,  a)  v’here  \>  is  the  total  number  of  customer!; 
served  by  station  I  and  o  la  the  total  number  of  customers  served  by  station 
II.  Each  time  a  customer  completes  service  at  either  station,  a  transition 
takes  place  to  state  (v  +  1,  o)  If  the  service  completion  occurred  at  station  I 
or  to  the  state  (v,a  +  1)  If  the  completion  occurred  at  station  11.  The  state 
(0,  0)  Is  the  starting  point  of  the  operation  of  the  queue  and  the  state 
(L,  L)  where  L  »  M  +  N  la  the  termination  point  of  the  operation,  since  all  of 
the  customers  have  been  served  at  this  time.  Each  realization  of  the  queueing 
process  Is  a  path  from  (0,  0)  to  (L,  L)  but  not  every  possible  path  Is  a 
possible  sequence  of  service  operations  by  our  queueing  system.  If  X(t)  and 
Y(t)  represent  the  number  of  customers  served  at  time  t  by  stations  I  and  II 
'..'respectively,  then  the  rules  of  operation  dictate  that 

(?)  X(t)  S  H  +  Y(t) 

Y(t)  <  N  +  X(t) 

These  conditions  determine  boundary  conditions  on  the  paths  from  (0,  0)  to 
(L,  L)  which  may  be  bona  fide  realizations  of  our  queueing  system.  Figure  1 
shows  the  boundaries  imposed  on  the  lattice  as  well  as  a  typical  path  represent¬ 
ing  a  realization  of  the  queueing  system.  The  realizations  are  constrained  to 
lie  within  an  Irregular  hexagon  H  determined  by  the  boundary  lines  X  ■  0, 

'  =»  0,  the  left  boundary  Y  =■  X  +  N,  the  upper  boundary  Y  =  M  +  N,  the  right 
boundary  X  =  M  +  N  and  the  lower  boundary  Y  *  X  -  M. 

’’Jlienever  a  cxistomer  finishes  service  at  station  I,  the  path  moves  to  thf. 
right  while  a  completion  at  station  II  causes  the  path  to  move  upward.  Whenever 
a  path  touches  either  the  left  boundary  or  the  lower  boundary,  the  path  must 
move  to  the  right  or  upward  respectively. 

4.  Computational  Methods  for  the  Total  Operation  Time  Distribution.  As 
pointed  out  above,  each  path  from  the  origin  (0,  0)  to  the  terminal  point  (L,  L) , 
'.;ii;'.ch  lies  wholly  within  the  hexagon  H  represents  an  operation  of  the  queueing 
system.  Also  associated  with  each  path  is  a  probability  distribution  on  the 
length  of  time  of  operation  of  the  queue  given  that  the  particular  path  was 
traversed.  With  this  In  mind,  we  may  describe  the  state  of  the  system  by  means 
of  a  quintuple 
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(2) 


(h,  j,  V,  T,  t) 


vhere  h  and  J  are  the  x,  y  coordinates  of  the  lattice  point  the  system 
currently  occupies,  and  are  the  number  of  units  of  service  time  given  to 

the  customers  at  stations  I  and  II  respectively  and  t  Is  the  total  time  elapsed. 

We  next  define  the  various  first  passage  probabilities  to  a  lattice 
point.  Suppose  that  the  process  Is  initially  In  state  (hg,  j^)  with  service  to 
the  customer  in  station  I  just  beginning  and  the  customer  at  station  II  hcvlng 
units  of  service  time  remaining  to  completion  of  his  service.  Assume  that  the 
system  arrives  for  the  first  time  In  state  (h,  j)  and  that  all  paths  from 
(hQi  Jg)  to  (h,  J)  are  permissible;  i.e.,  be  within  the  hexagon  H.  "'e  define  the 
quantity 

(3)  ^0’  t  -  0,  1, 

to  be  the  probability  that  t  units  of  time  are  required  for  a  passage  from 
(hg.  jg)  to  (h,  j)  where  the  passage  is  from  the  right  with  v  remaining  units 
of  service  required  for  the  customer  at  station  I  initially  and  t  tinits  of  time 
remaining  for  the  customer  at  station  II  when  the  passage  to  state  (h,  j) 
first  occurs.  Five  additional  probabilities  similar  to  \fi,  are  needed  the 
quantity 

U)  ^0* 

Is  similar  to  except  that  represents  the  number  of  remaining  units  of 
service  time  required  for  the  customer  at  station  II  rather  than  station  I 
when  the  state  of 'the  system  Is  (hg,  jg).  The  quantities  and  'l>2  are 
probabilities  of  transition  Into  a  state  where  the  transition  is  to  the  right 
into  the  final  state  (h,  j)  rather  than  from  below.  That  Is,  a  customer  at 
station  I  has  just  completed  his  service  while  the  customer  at  station  II  has  t  ' 
units  of  service  time  remaining. 
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Analogous  to  the  and  we  define 


1q»  h,  j,  V,  T,  t) 

to  be  the  probability  that  the  system  makes  a  transition  from  state  (h^,  j^) 
fo  the  state  (h,  j)  from  below  in  t  units  of  time  with  T units  of  service  time 
retoalning  for  the  customer  at  station  I  and  initially  v  units  of  time  remained  for 
for  the  customer  at  station  II  in  the  initial  state  (h^,  ig)*  Also  the  quantity 

^0* 

in  the  probability  of  transition  from  (h^,  j^)  to  (h,  j)  in  t  units  of  time  with 
T  service  units  remaining  for  the  customer  at  station  I  and  initially  v  units  of 
service  tine  remained  for  the  customer  at  station  I  rather  than  station  II. 

In  addition  to  the  transitions  from  the  right  (h  -  1,  j)-+-(h,j)  and 
the  transition  from  below  (h,  j  -  i)-*-(h,  j),  it  is  possible  to  have  a  diagonal 
transition  of  the  form  (h  -  1,  j  -  l)-^(h,j).  This  corresponds  to  two  customers 
one  at  each  station  both  finishing  service  simultaneously.  For  this  situation, 
we  need  two  additional  probabilities; 

(7)  ^0’ 

and 

(S)  e2(hQ,  Jq,  h,  j,  V,  t). 

Tliese  are  the  probabilities  that  a  diagonal  transition  from  state  (h^,  j^)  to 
rtate  (h,  j)  takes  place  in  t  units  of  time.  Additionally,  v  units  of  service 
time  remained  fro  the  customer  at  station  I  (0^^)  respectively  station  II  (62) 
in  the  initial  state  (h^,  Jq). 

Formulas  for  these  tr^sltion  probabilities  are  given  below. 

\ 

The  computation  of  the  total  service  time  distribution  proceeds  in 
four  separate  stages.  The  first  three  stages  consist  of  finding  the  probability 
distribution  on  the  time  of  first  passage  to  various  boundaries  in  the  hexagon 
n.  The  last  stage  consists  of  finding  the  passage  probability  to  the  terminal 
state  (L,  L). 


The  phase  I  boundary  consists  of  the  two  line  segment  joining  the 
points  (0.  N) ,  (M,  N)  and  (M,  0).  It  should  be  obvious  that  the  passage  from  the 
origin  (0,  0)  to  thi.,  boundary  occurs  with  probability  one.  The  probability 
distributions  at  lattice  points  along  this  boundary  are  given  by  the  quantities 

(9)  0.  M,  j.  0,  T,  t) 

where  j  -  0,  ....  N;  t  =  1,  . . . ,  L2 ;  and  t  -  0,  1,  T  where  T  is  the  maximum 

possible  total  service  time  required  to  reach  this  boundary.  These  correspond 
to  the  boundary  on  the  line  segment  from  (M,  0)  to  (M,  N)*  Along  the  segmen, 
from  (0,  N)  to  N)  the  probabilities  are  given  by 

(10)  <(>^(0,  0,  h,  N,  0,  T,  t) 

where  h  -  0,  ....  M.  t  =  1.  ....  and  t  -  0,  1,  ...  T.  The  comer  probability 
is  given  by 

(11)  0j^(O,  0,  M,  N,  0,  t). 

Note  that  since  v  -  0  in  each  case  the  subscript  1  could  be  replaced  by  a  2. 

The  phase  II  boundary  is  determined  by  the  line  segments  connecting 
the  points  (M.  L),  (M,  N)  and  (tl,  L).  The  computation  of  the  first  passage 
distribution  to  these  boundaries  is  considerably  more  complex  than  the 
probabilities  for  the  phase  I  boundary. 

The  lattice  points  along  the  phase  I  boundary,  constitute  a  partition 
of  the  probability  space.  Therefore,  by  the  law  of  total  probability,  the 
distribution  along  the  phase  II  boundary  may  be  obtained  by  computing  the 
distribution  given  that  the  process  started  at  a  point  along  the  phase  I 
boundary  and  summing  over  the  phase  I  boundary.  That  is 

(12)  P[(h,  j,  T,  t,  s)]  - 

Jp[(h,  j,  T,  t,  8)|(hQ,  Jq,  Tq,  tg,  8q) 1  P[(hQ,  Jq,  Tq,  t^,  Sq) 1 
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Figure  2  Boundaries  for  ’’rocessing 


where  (h,  j,  t  ,  t,  s')  is  a  state  along  the  phase  II  boundary,  (h,  j)  are  the 
lattice  points,  t  is  the  service  time  remaining  for  either  customer, 
s  =  1  or  2  denoting  the  customer  for  which  there  is  residual  service  time  and  t 
is  the  total  operational  time  up  to  this  point.  The  state  (h^,  j^,  Tq,  t^,  Sq) 
is  a  state  along  the  phase  I  boundary  with  the  quantities  h^,  Tq,  t^, 
having  the  same  meaning  as  described  for  the  phase  II  boundary.  The  summation 
is  taken  over  all  possible  values  of  (hg,  jg,  tg,  tg,  Sg) .  That  is,  for 
kg  =  M,  ,1^  -  0,  1,  ...,  N  -  1,  Tg  -  1,  2,  ...,  L^,  Bq  •  2  and  for  Jg  -  N, 
hg  *  0,  ...,  m  -  1,  Tg  ■  1,  2,  ...,  Ij^,  Sg  “  1.  Futther  details  of  this 
computation  are  given  in  the  section  on  computation. 

The  phase  III  computation  is  similar  in  strategy  to  the  phase  II 
strategy.  The  phase  III  boundary  consists  of  the  two  line  segments  connecting 
the  points  (M,  L) ,  (L,  L)  and  (L,  N).  The  strategy  is  to  regard  the  phase  II 
boundary  as  a  partition  of  the  probability  space  and  to  compute  the  conditional 
boundary  distribution  along  the  phase  III  boundary  given  that  the  process  passes 
through  the  particular  point  along  the  phase  II  boundary. 

The  phase  IV  computation  consists  of  finding  the  probability  distribu¬ 
tion  on  the  time  to  enter  the  state  (L,  L)  given  that  the  process  passes  through 
one  of  the  lattice  points  on  the  phase  III  boundary. 

Each  lattice  point  (h,  j)  along  the  phase  I  boundary  generates  a 
distribution  along  the  phase  II  boundary.  The  total  of  such  distributions  is 
M  +  N  +  1.  Each  lattice  point  of  this  distribution  generates  a  distribution 
along  the  phase  III  boundary.  Thus,  each  lattice  point  along  the  phase  I 
boundary  generates  (M  +  N  +  1) (M  +  H  -  1)  distributions  along  the  phase  III 
boundary.  These  distributions  are  joint  distributions  in  t  and  t.  Many  of 
these  distributions  have  a  total  probability  which  is  quite  small.  Thus, 
if  at  the  first  phase  of  processing,  some  distributions  are  eliminated  because 
the  passage  probability  at  that  point  is  small,  the  elimination  will  have  a 
negligible  effect  on  the  accuracy  of  the  total  operational  time  distribution, 
but  will  eliminate  the  computation  of  many  probability  distributions  which  will 
save  a  considerable  amount  of  computer  time. 
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The  large  number  of  distributions  involved  throughout  the  computation 
also  creates  a  need  for  a  substantial  amount  of  storage  for  these  distributions. 

T’o  solve  this  problem,  auxiliary  disk  storage  is  used  to  store  the  distributions 
along  the  particular  boundaries  for  each  phase  of  the  processing.  An  Indexing 
algorithm  was  developed  to  keep  tract  of  the  bookkeeping  Involved  In  this  process. 

A  check  on  computational  accuracy  and  correctness  la  made  at  each  stage 
of  the  processing  by  using  the  law  of  total  probability  on  each  of  the  boundaries 
for  phase  I,  phase  II  and  phase  III.  In  addition,  certain  Internal  auditing 
schemes  are  used  to  check  for  accuracy  and  correctness  of  the  algorithm. 

5.  Comp ut at  1  on al  Formulas  for  jj;,  ^  and  We  now  consider  the  problem  of 

computing  the  passage  probabilities  from  a  given  initial  lattice  point  (h^,  j^) 
to  another  lattice  point  (h»j).  It  is  assumed  that  the  process  remains  within 
the  square  whose  vertices  are  (h^,  Jq) .  (hp,  j),  (h,  j^)  and  (h,  j),  but  may 
otherwise  follow  any  path  from  (hp,  jp)  to  (h,  j).  We  consider  the  case  of  ’1'^^ 
first.  THe  passage  into  the  state  (h,  j)  is  from  the  fl^t.  In  addition,  v 
units  of  time  remained  to  service  the  customer  at  station  I.  A  transition  of 
this  type  takes  place  In  t  units  of  time  if  and  only  if  j  -  jp  customers  have  been 
served  at  station  II  and  service  has  begun  for  an  additional  customer  who  requires 
T  additional  units  of  service  to  complete  his  processing.  Simultaneously,  the 
customer  at  station  I  must  complete  his  v  units  of  service  time  and  an  additional 
h  -  hp  -  1  customers  must  be  served  with  the  last  customer  completing  his  service 
exactly  at  time  t. 

The  probability  that  v  customers  complete  their  service  at  station  II 
is  given  by  the  v-fold  convolution  of  the  station  II  service  time  distribution 

(13)  Y*''C) 

The  probability  that  the  customer  currently  in  service  requires  t  additional  units 
of  service  is  given  by 

(14)  Y(Tq  +  p 
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where  Tq  is  the  nunber  of  units  of  service  time  already  completed.  Therefore, 
the  probability  that  after  exactly  t  time  units,  the  station  II  condition  is 
satisfied  is 

Li-T 

(15)  Pjj  -  I  Y(Tq  +  T)  -  Tq) 

Tq-O 

Sindlarly,  we  may  derive  the  probability  Pj  that  the  station  I  conditions  are 
satisfied  is  the  probability  that  h  -  hp  -  1  customers  are  served  in  t  -  v  units 
of  time.  This  probability  la  given  by 

(16)  p^  -  -  v). 

Since  the  service  times  at  the  stations  are  presumed  to  be  Independent,  the 
probability  that  t  units  of  time  are  required  for  the  transition  is  given  by 

(17)  ’<'l(^0*  ^0’  “  ^(*^)  “  Pj  P£i  ■ 

L  — T 

X*(h'ho-l)(t  -  V)  -  Tq) 

To-0 

This  transition  probability  is  computed  for  each  value  of  t  ■  0,  1,  T 

where  T  is  the  largest  possible  value  of  time  during. which  a  transition  can 
take  place.  An  upper  bound  on  T  Is  given  by: 

(18)  T  -  min[1.2(J  -  j^),  v  +  Lj^(h  -  h^  -  1)]. 

In  the  actual  computational  procedure,  T  takes  on  a  considerably  smaller  value 
because  the  upper  and  lower  tails  of  the  probability  distributions  are  trimmed 
which  means  that  beyond  a  certain  point  In  the  distribution,  the  probabilities  are 
treated  as  though  they  were  zero.  Points  In  the  tails  of  the  distribution  which 
have  zero  probability  are  not  retained  in  the  storage  arrays.  This  greatly 
reduces  the  computation  time. 

By  arguments  analogous  to  the  above,  we  may  determine  the  formulas  for 
each  of  the  quantities  described  In  section  5.  We  omit  details. 
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(19) 


^0* 


r*(.h-h^) 


V" 


O^t)  I  Y(t  +  T)Y*^J'Jo"^\t  -  V  -  O 
^  =.n  ^  U 


To=>0 


(20) 


*^1^^0*  ^0*  '^* 


Y*CJ-jo)(t)  J  X(Tg  +T)X*^*'“V^^(t  -  V  -  Tq) 


Tq-O 


(21) 


^  2  ^^0  *  ^0*  ^  ^  ^  " 

Y*0-jo-l)(t  -  v)  ^  X(Tq  +  T)X*^^"^o\t  -  T.) 
Tq-O 


(22) 


'^0* 


Y*(J-Jo\t)X*^^‘V^^(t  -  v) 


(23)  <<’2(^0’  ^0’  '^’  '^’  " 

Y*(J~Jo~^\t  -  v)X*^^"^0^  (t) 

In  order  to  make  the  above  formulas  valid  for  all  J  ^  Jq  and  h  2  Hq,  we  use  the 
convention  that  for  any  probability  distribution  F,  F  (•)  is  the  probability 
distribution  which  assigns  probability  1  to  the  point  0. 

6 .  Computational  Procedures  for  Rectangular  Boundary  Exits.  The  phase  I 
processing  has  been  described  in  Section  4.  It  amotints  to  evaluating  the  formulas 
for  i|),  (^,  and  0  with  v  ■  0.  The  processing  algorithm  is  organized  so  that  the 
processing  begins  with  those  lattice  points  which  are  likely  to  have  high 
probability.  This  point  is  estimated  by  beginning  with  the  one  closest  to  the 
intersection  of  the  phase  I  boundary  with  the  line  through  the  origin  whose  slope 
is 
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(24) 


m 


4 

Zi^i 

1=1 

h 

1 1  \ 

i-1 


i.Q.,  the  ratio  of  the  means  of  the  Y  and  X  service  time  distributions. 


In  later  stages  of  processing,  it  Is  necessary  to  consider  an  evolution 
of  the  process  similar  to  the  phase  I  process.  For  example,  during  phase  III 
processing,  one  computes  the  conditional  distribution  along  the  phase  III 
boundary  given  that  the  process  passed  through  a  fixed  point  along  the  phase 
II  boundary,  ^^en  the  process  evolved  along  this  path  through  the  phase  II 
boxmdary,  one  of  the  customers  had  just  completed  his  service  while  the  other  had 
residual  service  time  remaining.  The  probability  distribution  along  the  phase  III 
boundaries  are  given  by  the  formulas 


(25)  P((h,  j,  T,  t)]  -  jg.  h,  j,  V,  T,t)P((hQ,  Jq,  V,  t,  1)) 

where  (h,  J)  Is  a  lattice  point  on  the  phase  III  boundary  and  is  the  residual 
service  time  for  ciistomers  at  station  I.  The  summation  Is  taken  over  all  lattice 
points  (hg,  jg)  along  the  phase  II  boundary,  ranges  over  Its  range 
(1,  ...,  If  1  *■  1  and  1,  . . . ,  L2  If  1  =  2)  and  1  =  1,  2.  This  formula  Is 
valid  for  probabilities  along  the  upper  boundary  connecting  the  points  (M,  L) 
and  (L,  L) . 


Along  the  right  boundary  connecting  the  points  (L,  N)  and  (L,  L)  the 
probability  Is  given  by 

(26)  P[(h,  j,  T,  t)]  -  Z'i’i(hg,  jg,  h,  j,  V,  T,  t)P[(hg,  jg,  V,  t,  1)] 

where  the  summation  extends  over  the  same  Indices  as  In  formula  (25). 


The  probability  distribution  at  the  comer  point  (L,  L)  Is  given  by 
(27)  P[(L,  L,  0,  t)]  -  Z^iChg,  jg,  L,  L,  V,  jg,  v,  t,  1)] 
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where  again  the  summation  extends  over  the  indices  defined  in  formula  (25) . 

These  computations  are  carried  out  by  a  subroutine  called  SQUARE. 

The  phase  IV  processing  for  total  operational  time  is  calculated  by  the 
same  procedure  as  for  the  phase  III  boundary,  except  that  it  is  a  degenerate  case* 
j  »  Jq  or  h  =  h^  depending  on  the  particular  segment  of  the  phase  IV  boundary. 

7.  Computational  Methods  for  Triangular  Boimdarles.  The  computational 
strategy  and  formulas  of  the  preceding  section  will  not  suffice  for  the  passage 
from  phase  I  to  phase  II  calculations  because  of  the  triangular  shaped  regloi/S 
in  which  the  process  must  evolve.  In  the  passage  from  the  phase  I  segment  defined 
by  the  points  (M,  0)  and  (M,  N)  to  the  phase  II  segment  defined  by  the  points 
(M,  N),  (I,  N) ,  the  condition  h  ^  J  +  M  must  be  observed  by  the  path.  VIhenever 
h  =  j  +  M,  idle  time  occurs  at  station  I  until  a  customer  finishes  at  station 
II.  This  customer  then  goes  to  station  I  for  service. 

The  phase  II  processing  can  be  reduced  to  a  sequence  of  rectangular 
boundary  problems  in  the  following  way.  Consider  the  lower  triangular  region. 

The  phase  I  distributions  are  computed  for  the  lattice  points  along  the  line 
segment  Joining  the  points  (M,  0)  and  M,  N).  At  the  point  (M,  0)  the  process 
must  make  a  transition  to  the  point  (M,  1);  no  other  transition  is  possible 
because  the  queue  at  station  I  is  now  empty  and  a  customer  at  station  II  must  be 

served.  From  any  point  along  the  line  segment  Joining  (M,  1)  and  (M,  N) ,  the 

process  must  pass  through  the  points  (H,  N) ,  (M  +  1,  N)  or  one  of  the  points 

along  the  segment  Joining  the  points  (M  +  1,  1)  and  (M  +1,  N).  The  exit 

distributions  at  these  points  are  calculated  using  the  same  methods  as  for  any 

rectangular  boundary.  The  rectangle  consists  of  the  segments  Joining  the  vertices 
(M,  1),  (M,  N),(M  +  1,  N),  (M  +1,  1).  This  rectangle  has  in  conmon  its  top 

side  (l.e.,  the  points  (M,  N)  and  (M  +  1,  N))  with  the  phase  II  boundary. 

The  processing  then  proceeds  with  the  rectangle  whose  vertices  are 
(M  +  1,  2),  (M  +  1,  N),  (M  +  2,  N),  (M  +  2,  2).  In  this  manner,  the  calculations 
are  made  step  by  step  with  the  phase  I  boundary  moving  in  the  x  direction  one 
step  at  a  time  until  the  lower  phase  I  boundary  is  mapped  onto  the  phase  II 
boundary. 
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The  upper  triangular  region  consisting  of  the  segments  connecting  the 
vertices  (0,  M) ,  (M,  N)  and  (M,  L)  is  processed  in  the  same  manner  as  the  lower 
triangular  region. 

8.  Tail  Trimming  and  Convolutions.  In  order  to  calculate  the  quantities 
4»,  4*,  and  §  the  convolutions  X  ■  0,  M  and  k  -  0,  ...,  N  are 

needed.  A  problem  arises  in  that  if  K  and  N  are  large,  a  large  amount  of  storage 
and  computational  time  la  needed  because  the  distributions  concentrate  on  a 
large  number  of  points.  Some  savings  can  be  effected  with  negligible  loss  in 
accuracy,  by  trimming  the  tails  of  the  higher  order  convolutions.  Even  though 
the  original  service  time  distribution  nay  concentrate  heavily  in  the  tails,  the 
tails  of  the  higher  order  convolutions  will  have  increasingly  smaller  probability 
because  of  the  central  limit  theorem  effect. 


An  adaptive  trimming  algorithm  vras  developed  to  permit  trimming  both 
the  upper  and  lower  tails  of  the  higher  order  convolutions.  Suppose  that  the 
probability  distribution  is  given  by 


(28) 


p^,  t  -  0,  1,  . . . ,  N 


Then  two  indices  t/>  and  t  are  chosen 


(29) 


corresponding  to  and  such  that 

<  Z  Pt 

t=0 


and 

N  N 

(30)  ^  p^  s  <  Jp^ 

The  distribution  is  then  approximated  by  the  new  distribution 


(31) 


Pt’ 


t  . 

u 
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The  tails  are  ignored.  The  effect  of  this  trimming  is  to  produce  a  defective 
distribution]  i.e.,  one  for  which  the  total  probability  is  less  than  one.  The 
total  probability  does,  however,  satisfy  the  inequality. 


(32) 


u 


I  Pt  ^  1  - 


'I 


e 

u 


1  -  e 


The  effect  on  the  mean  and  variance  of  this  trimming  process  can  easily  be 
estimated  analytically. 

l£t  y  be  the  mean  and  Uj.  be  the  mean  of  the  trimmed  distributions. 


t 


1 


(33) 


I  +  I  ip 

i=0  l=t 


1 


S  e»  t^  +  e  N  ^  Ne 

-C  -C  XX 


u 


similarly,  let  a  be  the  variance  and  Cj,  be  the  trimmed  variance.  Then 


'"I  N 

(33)  ^  I  i^p  +  I  i^p  -  (v  -  y.j.)(y  +  yj  s 

t=0  ^  t-t  ^  ^  ^ 

u 

2  2  2  2  2 
N  +  N  -  N  e  5  N  e 

The  effects  of  using  a  trimmed  distribution  on  convolutions  is  to  spread  the 
truncation  error  over  all  of  the  probabilities  rather  than  to  concentrate  it  at 
a  single  point.  Consider  a  distribution  F(*)  and  its  trimmed  version  F.j,(*)* 

Let  C  be  any  other  distribution.  Let  E(*)  ■  F(-)  -  F^(*)»  Then,  the  total  mass 
of  the  defective  distribution  F^(*)  la  |  1-e.  Consider  the  convolution 
F^*G.  The  total  mass  of  F^*G  is  I’^Clx*^),  the  product  of  the  total  mass  for  the 
two  distributions,  as  may  easily  be  shoi-m.  The  error  term  is  given  by 

(34)  F*G  -  F.j,*G  -  E*G 
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If  the  error  term  E 


Is  split  into  two  parts  the  part  trinmed  from  the  upper 
tail  and  the  part  trimmed  from  the  lower  tail,  then  we  may  consider  the  support  of 
E*G.  Suppose  for  simplicity  that  the  lower  tall  is  the  only  one  trimmed.  Then 
if  G  concentrates  on  the  points  0,  1,  ....  T„,  the  distribution  of  E*G 
concentrates  on  the  points  0,  1,  ^  maximum  probability  associated 

with  a  point  is  bounded  by  the  maximum  probability  at  any  point  in  the  G 
distribution  times  That  this  is  so  may  be  easily  seen  by  the  convolution 
formula 


H  H 

(35)  r_  -  t  q,)  I  P.  =  e  'T^ax  o, 

"  1=0  ^  ^  1=0 


where  q  is  the  value  of  the  probability  density  of  the  G  distribution  at  the 

ooint  y  and  r  is  the  probability  density  of  E*G  at  the  point  n. 

“■  n 

We  now  consider  the  savings  in  storage  and  the  loss  in  accuracy  due  to 
trimming.  A  specific  example  is  given.  For  this  example,  the  distribution  X 
is  uniform  on  10  points  and  vje  consider  the  convolutions  X  for  j  =  0,  1,  ...»  30. 
In  general,  the  storage  required  for  such  a  distribution  is  given  by  the  formula 

(36)  S  »  3  +  (L,  -  l)n(n  +  1) 

n  X  ~ 

where  L  is  the  number  of  points  in  the  distribution  for  X  and  n  is  the  number 
of  convolutions.  For  our  example  ■  10,  n  =  30,  and  *  4188. 

The  results  are  presented  in  Table  I.  Six  different  trimming  values  c 
are  given  along  with  the  total  number  of  support  points  for  all  thirty  distribu¬ 
tions  as  well  as  the  number  of  support  points  for  the  thirtieth  convolution  X 

*30 

The  mean  of  the  trimmed  version  of  X  as  well  as  the  total  probability  mass 
is  given.  The  trimming  procedure  used  was  to  trim  tails  after  each  convolution 
e  from  each  tail  of  the  distribution.  Thus,  for  the  fourth  distribution,  the 
following  was  formed: 


(((Xi.*X)^*X)^*X)^*X, 


where  the  subscript  T  denotes  trimming. 
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TABLE  1 

Effects  of  Trimming 


?riin  Value 

e 

Support  Points 
*30 

Total  X  ^ 

Mean 

Total  Mass 

0 

4188 

275 

165 . 0000 

1.00000000 

lO"^ 

3282 

175 

165.0000 

.99999997 

10-s 

3122 

163 

165.0000 

.99999969 

io“^ 

2939 

153 

165.0003 

.99999666 

10"^ 

2734 

139 

165.0022 

.99996895 

10-5 

2468 

123 

165.0218 

.99964830 

10-"^ 

2132 

103 

165.2055 

.99602595 
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9.  Bookkeeping  algorithm  for  Disk  Storage.  Bach  of  the  four  phases 
of  the  computation  results  in  the  computation  of  the  joint  distribution  of 
total  operational  time  and  residual  service  time  for  the  customer  being 
served  at  the  time  the  system  reaches  a  particular  lattice  point  on  the 
boundary.  These  distributions  are  stored  in  a  compressed  form  in  an  array. 
However,  because  of  the  large  amount  of  space  required  for  the  distributions, 
they  are  written  onto  random  access  disk  files  as  soon  as  the  distribution 
is  generated.  Only  one  lattice  point  at  a  time  is  processed  and  core  space 
for  the  distribution  at  only  one  point  is  provided.  As  the  processing 
proceeds  at  later  stages  e.g.  from  phase  I  to  phase  II,  it  is  necessary  to 
retrieve  distributions  from  various  points  along  the  boundary.  This 
necessitates  using  random  access  files  rather  than  sequential  access  files. 


The  standard  FORTRAN  random  access  procedures  are  utilized.  Hov/ever, 
this  access  method  involves  accessing  records  by  record  number.  The  boundary 
for  the  various  phases  of  the  processing  are  determined  by  lattice  points 
which  are  identified  by  pairs  (h,j).  A  translation  algorithm  was  devised 
for  translating  these  pairs  into  record  numbers.  The  algorithm  is  further 
complicated  by  the  fact  that  processing  along  any  boundary  is  not  necessarily 
in  a  logical  order  along  the  boundary.  Instead,  the  order  of  processing 
is  baaed  on  probability  values  in  order  to  take  advantage  of  tail  trimming. 
Another  complicating  factor  is  that  the  coordinate  pairs  for  each  boundary 
assume  different  values.  For  Instance,  along  the  phase  I  boundary,  the 
values  of  h  and  j  lie  in  the  intervals  OShSM.OSj^N.  On  the 
phase  II  boundary,  MShsM+N,  NSj  sM+N. 

Four  separate  files  are  utilized  for  the  processing.  Unit  1  is 
used  for  the  phase  I  boundary;  Unit  2  is  used  for  idle  time  information  along 
the  upper  triangular  boundary  and  the  lower  triangular  boundary;  Unit  3  is 
used  for  the  phase  II  boundary.  During  the  computation  of  the  phase  II 
boundary.  Unit  4  is  used  for  scratch  storage.  Unit  4  is  also  used  for  the 
phase  III  boundary. 

Corresponding  to  each  Unit,  an  array  is  established  with  each  position 
in  the  array  corresponding  to  a  lattice  point  along  the  boundary.  Initially, 
the  array  is  cleared  to  zeros  and  afterv/ard,  the  entry  in  the  array  is  the 
record  number  of  the  record  which  contains  the  distribution  for  that 
particular  lattice  point. 


Ihe  particular  records  for  a  lattice  point  are  addressed  by  a  pair 
(j,k)  where  j  =  1  or  2  representinQ  a  particular  segment  of  the  boundary 
and  k  is  either  the  x  or  y  coordinate  of  the  lattice  point  on  the 

boundary.  Thus,  for  phase  I,  the  points  along  the  segment  (h,N)  h  =  0,1,..., M 

are  addressed  by  the  pair  (l,h)  while  the  segment  (M,j),  j  =  0,1,..., N  is 
addressed  by  (2,j).  Similar  rules  are  used  for  the  other  units.  Table  II 
gives  the  array  index  for  the  lattice  points  along  the  phase  I  boundary, 

the  phase  II  boundary  and  the  idle  time  boundary.  The  example  given  has 

five  customers  initially  in  each  quene. 

Table  II  Indexing  System 


;  Unit 

1 _ 

1 

i  Unit 

1 

2  1 

Unit 

3 

h 

J 

index 

I 

1 

h 

i 

indek 

h 

j 

index 

5 

0 

1 

5 

0 

1 

5 

5 

6 

5 

1 

2 

' 

6 

5 

2 

6 

5 

5 

5 

2 

3 

7 

5 

3 

7 

5 

4 

5 

3 

4 

8 

5 

4 

8 

5 

3 

5 

4 

5 

9 

5 

5 

9 

5 

2 

5 

5 

6 

10 

5 

6 

10 

5 

1 

4 

5 

7 

0 

5 

7 

5 

6 

7 

3 

5 

8 

1 

6 

8 

5 

7 

8 

2 

5 

9 

2 

7 

9 

5 

8 

9 

1 

5 

10 

3 

8 

10 

5 

9 

10 

0 

5 

11 

1_ 

4 

9 

11 

5 

10 

11 

5 

10 

12 

10.  General  Results.  A  set  of  examples  was  chosen  to  illustrate  and  test 
the  computational  ideas  developed  here.  We  describe  briefly  the  characteristics 
of  each  example.  We  use  the  notation  (M,  N,  Lj^,  L^)  to  denote  the  particular 
example  where  M  is  the  number  of  customers  initially  at  station  I  and  Lj^ 
is  the  maximum  number  of  service  time  units  a  station  I  customer  will  need. 

At  station  II,  there  are  N  customers  with  a  maximum  number  of  service  units 
being 


The  following  examples  were  run. 

A.  (5, 5, 2,1).  The  service  time  at  station  II  required  one  unit  of 
service  with  probability  one.  The  other  station  had  a  uniform  service 
time  distribution. 

B.  (5, 5, 2, 2).  The  service  time  distribution  at  each  of  the  stations  vas 
uniform;  l.e.  with  probability  *s  one  unit  was  required  and  with 
probability  h  two  units  were  needed. 
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C.  (5, 5, 2, 3).  The  service  tine  distribution  at  station  I  is  unifonn; 

while  at  station  II,  =  .5,  =  .25,  =  .25 

D.  (5, 5, 3, 2).  This  example  is  the  same  as  example  3. 

E.  (5,10,2,2).  The  service  tine  distribution  is  uniform  at  each  station. 

F.  (10,5,2,2).  Again  the  service  time  distribution  is  uniform  at  each 
station. 

G.  (10,10,2,2).  The  service  time  distribution  is  uniform  at  each  station. 

H.  (15,15,4,4).  The  service  time  distribution  is  uniform  at  each  station. 

Example  H  is  the  largest  complete  computation  which  was  made.  Phase  I 
computations  were  made  for  an  example  with  30  customers  at  each  station. 

The  main  purpose  of  these  examples  was  to  serve  as  a  check  on  the  computations, 
test  the  program  for  correct  processing  of  asymmetric  examples  and  estimate 
processing  time. 

Table  111  lists  the  total  operational  time  distribution  for  the 
(5, 5, 2, I)  example.  The  distribution  is  not  binomial,  because  there  is  idle 
time  at  station  II  during  the  phase  II  processing. 

Table  III 

The  distribution  of  total  operational  time  for  the  (5, 5, 2,1)  system. 


■'  Time 

Probability  | 

Time 

Probability 

10 

.00097656 

16 

.17578125 

11 

.00976562 

17 

.10742187 

12 

.05371094 

18 

.04394531 

13 

.14648437 

19 

.00976562 

14 

.22460937 

20 

.00097656 

15 

i 

V.  .  ..  f 

.22656250 

The  complete  results  for  the  (15,15,4,4)  example  are  given  for  the 
un trimmed  case  and  the  trimmed  case.  Table  IV  gives  the  exit  probabilities 
at  the  phase  I  boundary.  Three  separate  calculations  are  given  to  show  the 
effects  of  trimming  on  these  probabilities.  In  the  untriramed  case,  no 
trimming  is  used,  while  in  the  second  case,  the  probabilities  are  trimmed 
at  e  =■  10  and  in  the  third  case,  the  probabilities  were  trimmed  at  10~®. 


j 
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Position 

.  H  K  ; 

_ _ 1 

Untrinir.-.^.  - 

PROBABILITY 

1  : 

£  =  iO"^ 

PROBABILITY  , 

r  =  10" 5 

PROBABILITY 

1 

15 

0 

i 

0.0 

\ 

i 

15 

1 

0.0 

1 

1 

15 

2 

0.0 

15 

3 

0.00000000 

15 

4 

0.00000001 

0.00000001 

15 

5 

0.00000091 

0.00000090 

0.00000074 

15 

6 

0.00002313 

0.00002313 

0.00002288  i 

15 

7 

0.00026704 

0.00026704 

0.00026680 

15 

8 

0.00175876 

0.00175875 

0.00175851 

15 

9 

0.00753842 

0.00753842 

0.00753822 

15 

10 

0.02287019 

0.02287019 

0.02287004 

15 

11 

0.05202916 

0.05202916 

0.05202888 

15 

12 

0.09255192 

0.09255192 

0.09255180 

15 

1 

13 

0.13235647 

0.13285646 

[  0.13285630 

15 

14 

0.15771523 

0.15771523 

0.15771513 

15 

15 

0.06477752 

0.06477752 

j  0.06477745 

0 

15 

0.0  1 

1 

0.00000001 

0.00000074 

1 

15 

0.0 

0.00000090 

0.00002288 

2 

15 

0.0 

0.00002313 

0.00026680 

3 

15 

0.00000000 

0.00026704 

0.00175851 

4 

15 

0.00000001 

0.00175875 

0.00753822 

5 

15 

0.00000091 

0.00753842 

0.02287004 

6 

15 

0.00002313 

0.02287019 

0.05202888 

7 

15 

0.00026704 

0.05202916 

0.09255180 

8 

15 

0.00175876 

0.09255192 

0.13285630 

9 

15 

0.00753842 

0.13285646 

0.15771513 

10 

15 

0.02287019 

0.15771523 

11 

15 

0.05202916 

! 

1 

12 

15  i 

0.09255192 

1 

1 

13 

15 

0.13285647 

i 

14 

1 

15 

0.15771523 

i~-  ----- 

1 

Table  IV 

Phase  I  exit  probabilities 
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Trimming  is  understood  to  mean  that  in  distributions,  tail  probabilities 
which  are  below  e  (=  10~®  and  10  ®  respectively)  are  truncated  to  zero  and 
these  points  are  deleted  from  the  distribution.  Moreover,  any  conditional 
distribution  at  any  of  the  lattice  points  in  the  hexagon  H  whose  total 
probability  is  less  than  e  are  completely  ignored. 

Table  Va  gives  the  total  operational  time  distribution  for  the 
untrimmed  and  each  of  the  two  trimmed  examples.  Table  Vb  gives  a  summary  of 
the  computational  time  for  each  phase  of  the  computation.  Table  Vc  gives 
a  list  of  the  amount  of  probability  lost  at  each  stage  due  to  trimming. 

Table  VI  gives  the  probability  distribution  at  the  phase  III  boundary.  These 
are  closely  related  to  the  idle  time  distributions. 


Table  Va 

Operational  Time  Distribution 


Time 

e  =  0 

Probability 

e  =  10-® 
Probability 

e  “  10-6 

Probability 

55 

0.0 

0.0 

0.0 

56 

0.000001 

0.000001 

0.000001 

57 

0.000003 

0.000002 

0.000002 

58 

0.000007 

0.000006 

0.000006 

59 

0.000019 

0.000016 

0.000016 

60 

0.000046 

0.000041 

0.000041 

61 

0.000108 

0.000096 

0.000097 

62 

0.000239 

0.000213 

0.000215 

63 

0.000500 

0.000448 

0.000453 

64 

0.000990 

0.000893 

0.000903 

65 

0.001859 

0.001687 

0.001705 

66 

0.003312 

0.003024 

0.003055 

67 

0.005608 

0.005148 

0.005202 

68 

0.009029 

0.008334 

0.008421 

69 

0.013838 

0.012837 

0.012971 

70 

0.020202 

0.018832 

0.019030 

71 

0.028120 

0.026334 

0.026614 

72 

0.037350 

0.035130 

0.035507 

73 

0.047378 

0.044741 

0.045228 

74 

0.057442 

0.054445 

0.055054 

75 

0.066623 

0.063358 

0. 064085 

76 

0.073981 

0.070563 

0.071400 
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77 

0.078724 

0,075275 

0.076204 

78 

0.080348 

0.076986 

0.077980 

79 

0.078726 

0.075549 

0.076576 

80 

0.074123 

0.071203 

0.072228 

81 

0.067126 

0.064509 

0.065499 

82 

0.058526 

0.056233 

0.057158 

83 

0.049174 

0.047207 

0.048043 

84 

0.039852 

0.038200 

0.038931 

85 

0.031180 

0.029820 

0.030439 

86 

0.023569 

0.022475 

0.022982 

87 

0.017225 

0.016366 

0.016768 

88 

0.012178 

0.011521 

0.011829 

89 

0.008333 

0.007843 

0.008072 

90 

0.005520 

0.005166 

0.005331 

91 

0.003541 

0.003292 

0.003407 

92 

0.002199 

0.002030 

0.002107 

93 

0.001322 

0.001211 

0.001261 

94 

0.000769 

0.000699 

0.000730 

95 

0.000433 

0.000390 

0.000409 

96 

0.000235 

0.000210 

0.000221 

97 

0.000123 

0.000109 

0.000115 

98 

0.000063 

0.000055 

0.000058 

99 

0.000031 

0.000026 

0.000028 

100 

0.000014 

0.000012 

0.000013 

101 

0.000006 

0.000005 

0.000006 

102 

0.000003 

0.000002 

0.000003 

103 

0,000001 

0.000001 

0.000001 

Table  Vb 

computation  time 


Phase 

e  =  0 

e  =  10"® 

e  -  10"® 

I 

1.016 

.958 

.923 

11 

38.341 

27.723 

21.526 

III 

157.054 

110.413 

96.315 

IV 

1.376 

1.037 

.926 

Total 

197.787 

140.131 

119.690 

Table  Vc 


Trimming  Effects 


Probability  lost 

Phase 

£  »  10"® 

e  =  10"® 

I 

0 

3  X  10~® 

4  X  10"® 

11 

0 

3  X  10"2 

1  X  10"^ 

III 

0 

2  X  10"2 

2  X  10-2 

IV 

0 

0 

0 

The  total  number  of  points  on  which  the  final  operational  time 
distribution  resides  is  95  for  e ■  0,  68  for  e *  10  ®  and  62  for  e  ■  10 
A  strange  pattern  exists  between  the  e  values  and  the  probability  mass 
lost  due  to  trimming  The  case  £ ■  10  ®  seems  to  have  lost  more  probability 
than  the  case  e«  lO”®.  Also,  the  total  probability  lost,  particularly  in 
the  later  stages  o£  the  computation,  is  relatively  large.  A  smaller  ^ 
value  in  these  late  stages  of  the  computations  would  probably  remedy  the 
situation. 
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Table  VI 

Phase  III  boundary  probabilities 


H 

K 

e  »  0 

e  =  10"® 

e  -  10"® 

30 

15 

0.00000047 

30 

16 

0.00000353 

0.00000003 

0.00000093 

30 

17 

0.00001917 

0.00000217 

0.00001085 

30 

18 

0.00008556 

0.00002586 

0.00005624 

30 

19 

0.00031799 

0.00015313 

0.00022866 

30 

20 

0.00100092 

0.00060582 

0.00076979 

30 

21 

0.00271011 

0.00184516 

0.00217903 

30 

22 

0.00640707 

0.00472150 

0.00536664 

30 

23 

0.01341559 

0.01062907 

0.01178075 

30 

24 

0.02516235 

0.02135680 

0.02314250 

30 

25 

0.04252185 

0.03824821 

0.04054072 

30 

26 

0.06477119 

0.06079866 

0.06320016 

30 

27 

0.08875024 

0.08566056 

0.08772284 

30 

28 

0.10925278 

0.10721851 

0.10868396 

30 

29 

0.12099172 

0.11984512 

0.12072813 

30 

30 

0.04917891 

0.02668752 

0.02668738 

15 

30 

0.00000047 

0.00000087 

0.00000275 

16 

30 

0.00000353 

0.00001266 

0.00001852 

17 

30 

0.00001917 

0.00007661 

0.00008490 

18 

30 

0.00008556 

0.00031027 

0.00031719 

19 

30 

0.00031799 

0.00099642 

0.00100009 

20 

30 

0.00100092 

0.00270824 

0.00270923 

21 

30 

0,00271011 

0.00640649 

0.00640630 

22 

30 

0.00640707 

0.01341545 

0.01341467 

23 

30 

0.01341559 

0.02516230 

0.02516173 

24 

30 

0,02516235 

0.04252140 

0.04252101 

0.06477032 

0.08874692 

0.10923414 

0.12091419 


25 

26 

27 

28 
29 


30 

30 

30 

30 


0.04252185 
0.06477119 
0.08875024 
!  0.10925278 


0.06476609 

0.08871662 

0.10910787 

0.12054744 


30  ;  0.12099172 
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11.  Uaiting  time  for  the  ith  customer.  The  problem  of  finding  the  waiting 
time  distribution  for  the  ith  customer  initially  at  station  I  (or  station  II) 
is  easily  solved  using  the  method  and  programs  described  in  the  previous 
sections.  Suppose  that  we  wish  to  find  the  waiting  time  distribution  of  the 
customer  who  is  initially  customer  i  at  station  I.  This  customer  i  must 

be  served  by  both  stations  I  and  II.  All  of  the  customers  in  line  behind 
customer  i  will  necessarily  be  served  after  he  is.  As  far  as  customer  1  is 
concerned,  these  customers  have  no  relevance.  Thus  one  can  consider  a  reduced 
queueing  system  which  has  initially  i  customers  at  station  I  and  N  customers 
at  station  II. 

In  the  reduced  queueing  system,  customer  i  is  now  the  last  customer 
at  station  I.  His  total  waiting  time  is  that  required  for  him  to  clear 
station  II.  We  may  describe  the  waiting  time  problem  in  terras  of  the  random 
walk  representation  pictured  in  figure  1.  In  this  representation,  the  last 
customer  initially  at  station  I  clears  station  II  vjhen  the  sojourn  path 
beginning  at  the  origin  reaches  the  upper  boundary,  i.e.  the  segment  connecting 
the  lattice  points  (i,L)  to  (L,L).  At  this  time,  station  I  may  still  be 
serving,  but  station  II 's  service  is  complete. 

In  the  notation  of  the  preceding  sections,  a  formula  may  be  written 
for  the  waiting  time  distribution  for  the  ith  customer  as  follows. 

N  Lj 

(37)  p  -  I  I  P(i  +  k,  L,  T,  t),  t  -  0,1 . 

k*l  T=1 

where  P  is  the  probability  distribution  at  the  lattice  point  (i  +  k,  L) 
along  the  upper  boundary,  and  x  is  the  number  of  service  time  units  remaining 
for  the  customer  at  station  I. 

12.  Idle  time  distribution.  There  are  two  types  of  idle  time  at  each  of 
the  serving  stations.  At  station  I,  idle  time  occurs  at  the  lower  boundary 
in  figure  1  which  corresponds  to  the  line  at  station  I  being  empty,  but  there 
remain  customers  at  station  II  who  must  also  be  served  by  station  I.  Idle 
time  also  occurs  at  the  right  boundary  in  figure  1.  This  corresponds  to  the 
completion  of  service  at  station  I  for  all  customers  in  the  system,  but  there 
yet  remain  customers  at  station  I  who  require  service.  Similarly  at  station  II 
idle  time  of  each  type  occurs. 
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The  idle  time  distribution  of  the  second  type  is  relatively  easy  to 
obtain  once  the  probability  distribution  for  the  phase  III  boundary  is  obtained. 
This  distribution  is  given  for  station  I  by 

L  L2 

P(t)  -  I  I  y*^^  “  (t  -  t)  P  [j.  L,  t] 
j-N  T-1 

where  P  [j,  L,  t]  Is  the  probability  of  passage  to  the  lattice  point  (J,L) 
with  T  units  of  residual  service  time  at  station  II.  This  probability  is 
easily  obtained  from  the  phase  III  boundary  probabilities  by  summing  out  the 
time  parameter.  A  similar  formula  is  available  for  the  station  II  idle  time 
distribution. 

The  other  type  of  idle  time  distribution  v.hich  occurs  at  the  lower 
boundary  and  at  the  left  boundary.  This  probability  calculation  is  not  as 
readily  available  as  a  by  product  of  the  total  operational  time  distribution. 

The  programs  and  algorithms  could  be  modified  to  keep  tract  of  the  idle  time 
of  this  type  as  well  as  the  total  operation  time.  This  procedure  Xijould  amount 
to  keeping  tract  of  the  joint  distribution  of  the  idle  time  as  v;ell  as  the 
total  operational  time.  This  would  require  a  considerable  Increase  in  the 
amount  of  storage  needed.  Another  alternative  is  to  modify  the  algorithm 
to  keep  tract  only  of  the  idle  time  at  the  triangular  boundaries. 
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II.  Cor.parison  VJith  Sir.iulation  Methods 

A  simulation  study  was  made  in  order  to  assess  the  exact  computational 
algorithms  for  correctness.  An  additional  benefit  was  to  determine  how  much 
accuracy  in  the  operational  time  distribution  could  be  obtained  by  simulation 
methods.  In  this  section,  we  discuss  the  random  number  generation  methods 
and  testing,  the  simulation  algorithm,  and  a  collection  of  results. 

Random  number  generation.  Three  methods  for  generating  random  numbers 
uniformly  distributed  over  the  unit  interval  were  considered.  They  were: 

(a)  linear  congruential,  (b)  random  numbers  from  an  adapted  form  of  the 
Rand  Corporation's  million  random  digits,  and  (c)  a  mixture  of  the  two. 

The  linear  congruential  generator  used  was  the  usual  one  in  which  the 
recurrence  formula 

X  ,  .  =  aX  mod(b) 
n+1  n 

is  used  to  generate  the  random  sequence.  In  this  problem,  the  value  of  a 

-JO 

used  was  32771  and  the  value  of  b  vtas  2  . 

The  random  numbers  from  the  Rand  Corporation's  million  random  digits 
were  converted  to  twenty-four  bit  binary  Integers  and  used  as  random  numbers 
distributed  uniformly  on  the  unit  interval  by  a  suitable  positioning  of  the 
binary  point. 

One  of  the  criticisms  of  linear  congruential  generators  is  the 
possibility  of  serial  correlations  and  lack  of  randomness  in  runs.  The  purely 
random  numbers  suffer  from  the  relatively  small  number  of  them.  In  our  case 
the  total  number  of  such  nurjbera  wa»  about  136,000. 

A  method  suggested  for  solving  these  difficulties  is  to  mix  the  random 

numbers  with  pseudo  random  numbers  in  the  following  way.  Generate  a  sequence 

of  pseudo  random  numbers  ^X^}  as  described  above  as  x/ell  as  a  sequence  of 

random  numbers  {Y  }  which  have  the  property  that  Y  +  138,000  -  Y  ,  i.e.  use 
n  r  r  j  n  n 

the  first  138,000  random  numbers  and  after  that  reuse  the  original  sequence. 
The  mixed  random  number  Z  is  obtained  by  Z  •■X  9  Y  v/here  the  symbol  9 
means  that  Z^  is  obtained  using  the  exclusive  or  on  the  digits  of  the  binary 
representation  of  and  Y^  to  obtain  the  binary  digits  of  Z^.  The  advantage 
of  this  mixing  is  that  the  sequence  enjoys  the  randomness  properties  of 

the  random  numbers  and  also  the  advantages  of  the  pseudo  random  sequence. 
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For  each  of  the  nethods  of  random  number  generation  described  above, 
the  following  tests  were  made:  (a>  distribution  test,  (b)  serial  correlation 
test,  and  (c)  runs  test.  For  each  test,  10,000  random  numbers  were  generated. 
The  distribution  test  consisted  of  partitioning  the  unit  interval  into  1000 
Subintervals  of  equal  length,  obtaining  the  empirical  frequency  count  and 
applying  a  chi-square  goodness  of  fit  test.  The  serial  correlation  test 
consisted  of  computing  the  one  step  serial  correlation  for  the  pairs 
Y  =  1,2,...,  n  -  1.  The  runs  test  is  described  in  HI]  p.  60.  In  this 
particular  test,  we  counted  runs  of  length  1,2, 3, 4, 5  and  greater  than  or  equal 
to  6.  The  observed  number  of  runs  was  compared  with  the  expected  and  a  chi- 
square  value  was  calculated.  In  all  cases,  the  chi-square  values  were  close 
to  the  expected  values  as  well  as  the  overall  mean  and  variance.  The  results 
are  given  in  table  VII.  They  suggest  that  any  of  the  methods  of  random  number 
generation  are  acceptable. 


Table  VII 
Simulation  Results 

_ pseudo  random  random _ mixed 


overall  mean 

.5030 

.5008 

.4978 

overall  variance 

.0838 

.0831 

.0837 

serial  correlation 

.003644 

-.001877 

-.004844 

distribution  chi  square  (999df) 

1036.6 

1028.4 

1029.4 

run  length  chi  square  (6df) 

5.471 

7.136 

4.092 

Simulation  Algorithm.  A  relatively  straight  fon^ard  simulation  algorithm 
was  used.  For  a  given  queue  setup,  a  set  of  service  times  was  generated  for 
each  station,  according  to  the  service  time  distribution  for  that  station, 
and  for  each  of  the  M  +  N  customers  to  be  served  by  that  particular  station. 

The  successive  service  times  were  independent  from  customer  to  customer  as 
Well  as  from  statfon  to  station.  Once  the  service  times  v;ere  generated,  an 
algorithm  determined  the  overall  operational  time  for  the  system.  The  total 
service  time  was  not  merely  the  larger  of  the  sum  of  the  two  service  times 
for  each  station,  because  of  the  possibility  of  idle  time  at  one  station  while 
it  awaits  customers  from  the  other.  This  consideration  was  easily  accomodated 
by  the  algorithm.  We  omit  further  details. 

3-  Comparison  of  Numerical  Results  with  Simulation  Results.  We  now  list 
for  comparison  some  simulation  results  with  exact  numerical  computations.  A 
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chi  square  goodness  of  fit  test  was  made  to  compare  the  numerical  results  with 
the  simulation  results  for  the  total  operational  time  distribution  of  examples 
B  through  H  of  section  I.  10.  The  chi  square  values  are  given  in  table  VIII. 
In  table  IX,  the  total  operational  time  distribution  for  the  (5,  5,  2,  3) 
example  is  given  as  obtained  by  simulation.  The  number  of  simulations  was 
10,000  in  each  case. 

Table  VTII 

Comparison  of  Simulation  With  Exact  Calculations 

Empirical 

True  Distribution  Distribution  . 


Example 

Chi  square  df 

Mean 

Variance 

Mean 

Variance 

(5,  5, 

2, 

2) 

31.34 

9 

15.84 

1.67 

15.90 

1.74 

(5.  5, 

2, 

3) 

93.10 

15 

17.75 

4.82 

17.88 

5.34 

(5.  5, 

3, 

2) 

51.07 

15 

17.75 

4.82 

17.86 

5.09 

(10,  5, 

2, 

2) 

14.49 

11 

23.62 

2.66 

23.69 

2.57 

(10,10, 

2, 

2) 

24.92 

13 

31.20 

3.34 

31.27 

3.36 

(15,15, 

4, 

4) 

45.44 

33 

78.27 

24.79 

78.47 

25.62 

The  general  results  show  that  as  far  as  the  mean  and  variance  are  concerned, 
the  simulation  produces  close,  but  not  accurate  results.  However,  as  far  as 
the  overall  distribution  is  concerned,  the  good  ness  of  fit  is  poor,  especially 
when  there  are  a  small  number  of  customers  in  the  system.  An  analysis  of  the 
results  shows  that  the  poor  fit  is  in  the  tails  of  the  distribution.  An 
extreme  example  is  given  in  Table  IX  for  the  (5,  5,  2,  3)  example  which  gave  the 
poorest  fit  of  all  the  examples  run. 
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Table  IX 

Goodness  of  Fit  Results 


Operational 

Time 

Exact 

Probability 

Frequency 

Count 

Expected 

Frequency 

Contribution 
to  x2 

10 

0.00000095 

0 

0.0 

0.01 

11 

0.00006199 

0 

0.6 

0.62 

12 

0.00110626 

10 

11.1 

0.10 

13 

0.00869751 

105 

87.0 

3.74 

14 

0.03723145 

378 

372.3 

0.09 

15 

0.09591246 

923 

959.1 

1.36 

16 

0.16157508 

1517 

1615.8 

6.04 

17 

0.19060600 

1869 

1906.1 

0.72 

18 

0.17199892 

1621 

1720.0 

5.70 

19 

0.13040411 

1259 

1304.0 

1.56 

20 

0.08936733 

957 

893.7 

4.49 

21 

0.05579908 

614 

558.0 

5.62 

22 

0.03127172 

378 

312.7 

13.63 

23 

0.01547158 

227 

154.7 

33.77 

24 

0.00676496 

97 

67.6 

12.73 

25 

0.00257871 

35 

25.8 

3.29 

26 

0.00085335 

7 

8.5 

0.28 

27 

0.00023588 

2 

2.4 

0.05 

28 

0.00005317 

1 

0.5 

0.41 

29 

0.00000861 

0 

0.1 

0.09 

30 

0.00000090 

0 

0.0 

0.01 
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III.  General  Conclusions 

The  results  of  this  study  have  shown  that  it  is  possible  to  analyze 
complex  queueing  models  of  the  type  discussed  in  this  report  using  discrete 
time  methods  together  with  modem  computational  facilities.  Heretofore, 
analytic  methods  have  been  able  to  solve  only  the  simplest  cases  of  queues 
having  exponential  service  times.  The  analytic  methods  rely  heavily  on  the 
memory less  property  of  the  exponential  distribution.  Using  discrete  time 
methods,  proposed  by  Neuts  [10]  and  others,  the  queueing  problem  is  completely 
solvable  with  arbitrary  service  time  distributions.  The  results  reported  here 
show  that  queues  with  fifteen  customers  initially  at  each  station  and  a  service 
time  distribution,  a  very  small  amount  of  computer  time  is  needed.  With 
service  time  distributions  concentrating  on  a  larger  number  of  points,  say 
sixteen  points,  for  each  station,  the  computation  time  would  Increase  by  a 
factor  of  sixteen  or  less  depending  on  the  effects  of  trimming.  The  amount 
of  computer  time  needed  would  still  remain  within  reasonable  bounds  —  about 
thirty  minutes  based  on  extrapolations  from  out  results.  Still  larger 
problems  could  be  solved,  but  refinements  of  the  algorithms  would  be  needed. 

In  particular,  the  use  of  the  Fast  Fourier  Transform  algorithm  for  convolutions 
of  distributions  concentrating  on  a  large  number  of  points  would  considerably 
decrease  processing  time.  For  the  small  examples  studied  in  this  report,  the 
use  of  the  Fast  Fourier  Transform  would  not  be  efficient,  since  convolution 
via  this  transform  method  is  faster  than  the  direct  methods  only  when  the 
particular  distributions  Involved  concentrate  on  more  than  128  points. 

The  commonly  used  simulation  methods  for  solving  queueing  problems  of 
this  type  are  not  as  useful  for  obtaining  the  operational  time  probability 
distribution  as  their  common  use  would  Indicate.  Our  experience  indicates 
that  the  mean  and  variance  can  be  obtained  reasonably  well  if  one  is  satisfied 
with  two  digit  accuracy.  An  additional  digit  or  so  is  possible  by  increasing 
the  number  of  simulations.  However,  the  frequency  counts  do  not  give  accurate 
estimates  for  the  probability  density.  The  problem  is  more  serious  for  queues 
with  smaller  numbers  of  customers  rather  than  larger  ones.  The  chi  square  tests 
generally  fall  until  there  are  ten  customers  in  the  queue.  This  phenomena  la 
probably  due  to  a  tendency  for  the  operational  time  distribution  to  exhibit  a 
limiting  behavior  as  the  number  of  customers  at  a  station  Increases. 


-33- 


One  of  the  major  drawbacks  to  this  method  of  analysis  is  the  rather  large 
amount  of  complex  computer  programming  which  is  required.  The  overall  amount  of  t. 
time  involved  is  certainly  larger  than  that  required  for  a  simulation,  but 
probably  not  much  more  than  would  be  required  via  analytic  methods. 

A  further  aspect  of  this  problem  which  needs  to  be  analyzed  is  the 
approximation  procedure  which  would  be  required  to  model  an  exponential  service 
time  distribution  by  a  discrete  time  distribution,  say  a  truncated  geometric 
distribution. 
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